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ABSTRACT 

Many astrophysical processes involving magnetic fields and quasi-stationary processes are 
well described when assuming the fluid as a perfect conductor. For these systems, the ideal- 
magnetohydrodynamics (MHD) description captures the dynamics effectively and a number 
of well-tested techniques exist for its numerical solution. Yet, there are several astrophysical 
processes involving magnetic fields which are highly dynamical and for which resistive effects 
can play an important role. The numerical modeling of such non-ideal MHD flows is signifi- 
cantly more challenging as the resistivity is expected to change of several orders of magnitude 
across the flow and the equations are then either of hyperbolic -parabolic nature or hyperbolic 
with stiff terms. We here present a novel approach for the solution of these relativistic resistive 
MHD equations exploiting the properties of implicit-explicit (IMEX) Runge Kutta methods. 
By examining a number of tests we illustrate the accuracy of our approach under a variety of 
conditions and highlight its robustness when compared with alternative methods, such as the 
Strang-splitting. Most importantly, we show that our approach allows one to treat, within a 
unified framework, both those regions of the flow which are fluid-pressure dominated (such 
as in the interior of compact objects) and those which are instead magnetic -pressure domi- 
nated (such as in their magnetospheres). In view of this, the approach presented here could 
find a number of applications and serve as a first step towards a more realistic modeling of 
relativistic astrophysical plasmas. 
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1 INTRODUCTION 

A vast number of astronomical observations suggests that mag- 
netic fields play a crucial role in the dynamics of many phenonema 
of relativistic astrophyics, either on stellar scales, such as for pul- 
sars, magnetars, compact X-ray binaries, short and long/gamma- 
ray bursts (GRBs) and possibly for the collapse of massive stellar 
cores, but also on much larger scales, as it is the case for radio 
galaxies, quasars and active galactic nuclei (AGNs). A shared as- 
pect in all these phenomena is that the plasma is essentially electri- 
cally neutral and the frequency of collisions is much larger than the 
inverse of the typical timescale of the system. The MHD approxi- 
mation is then an excellent description of the global properties of 
these plasmas and has been employed with success over the several 
decades to describe the dynamics of such systems well in their non- 
linear regimes. Another important common aspect in these systems 
is that their flows are characterized by large magnetic Reynolds 
numbers TZm = LV/X = AnaLV / c 2 , where L and V are the 
typical sizes and velocities, respectively, while A is the magnetic 
diffusivity and a is the electrical conductivity. For a typical rela- 
tivistic compact object, TZyi 1 and, under these conditions, the 



magnetic field is essentially advected with the flow, being contin- 
uosly distorted and possibly amplified, but also essentially not de- 
caying. We note that these conditions are very different from those 
traditionally produced in the Earth's laboratories, where TZyi -C 1, 
and the resistive diffusion represents an important feature of the 
magnetic-field evolution. 

A particularly simple and yet useful limit of the MHD ap- 
proximation is that of the "ideal-MHD" limit. This is mathe- 
matically defined as the limit in which the electrical resistiv- 
ity 77 = 1/(7 vanishes or, equivalently, by an infinite electri- 
cal conductivity. It is within this framework that many multi- 
dimensional numerical codes have been developed over the last 
decade to study a number of phen omena in relativistic astrophysics 
and in fully nonlinear regimes i Komissarov 1999b; Koid e et all 
19991: lKomissarovll200ll ; [Koldoba et al . 2002; Ga mmie et al.l | 2003|: 



Del Zanna et all 120031; lArminos et all 120051: i Duez et all J2005 
Shibata & Sekiguchill2005t iNeilsen et al.ll2006l: Unton et all 20061 ; 
McKinnevI l2006al; iMignone & Bodol l2006r. iNoble et alj 120071: 
Giacomazzo & Rezzollall2007l ; iDel Zanna et alj|2007l ; iFarris etltU 
20081) . The ideal-MHD approximation is not only a convenient way 
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of writing and solving the equations of relativistic MHD, but it is 
also an excellent approximation for any process that takes place 
over a dynamical timescale. In the case of an old and "cold" neu- 
tron star, for example, the electrical and thermal transport prop- 
erties of the matter are mainly determined by the transport prop- 
erties of the electrons, which are the most important carriers of 
charge and heat. At temperatures above the crystallization tem- 
perature of the ions, the electrical (and thermal) conductivities are 
governed by electron scattering off ions and an app roximate ex- 
pression for the electrical conductivity is given by dLamblll99lh 
g « 10 24 (10 9 K/T) 2 0/10 14 gem" 3 ) 3/4 s" 1 , where T and 
p are the stellar temperature and mass densitjQ. Even for a mag- 
netic field that varies on a length-scale as small as L ~ 0.1R, 
(where R is the stellar radius) the magnetic diffusion timescale is 



Tdiff = <lirL a/c 



3 x lCTyr- 



Clearly, at these temperatures and densities, Ohmic diffu- 
sion will be neglible for any process taking place on a dynami- 
cal timescale for the star, i.e., < few s, and thus the conductiv- 
ity can be considered as essentially infinite. However, catastrophic 
events, such as the merger of two neutron stars, or of a neutron star 
with a black hole, can produce plasmas with regions at much larger 



temperatures (e.g., T 



10 1 



K) and much lower densities 



(e.g., p ~ 10 8-10 gem -3 ). In such regimes, all the transport prop- 
erties of the matter will be considerably modified and non-ideal 
effects, absent in perfect-fluid hydrodynamics (such as bulk vis- 
cosity) and ideal MHD (such as Ohmic diffusion on a much shorter 
timescale Tdis ~ 10 3 s) will need to be taken into account. Similar 
conditions are likely not limited to binary mergers but, for instance, 
be present also behind processes leading to long GRBs, thus ex- 
tending the range of phenomena for which resistive effects could be 
important. Note also that these non-ideal effects in hydrodynamics 
(MHD) are proportional not only to the viscosity (resistivity) of the 
plasma, but also to the second derivatives of the velocity (magnetic) 
fields. Hence, even in the presence of a small viscosity (resistivity), 
their contribution to the overall conservation of energy and momen- 
tum can be considerable if the velocity (magnetic) fields undergo 
very rapid spatial variations in the flow. A classical example of the 
importance of resistive MHD effects in plasmas with high but fi- 
nite conductivities is offered by current sheets. These phenomena 
are often observed in the solar activity and are responsible for the 
reconnection of magnetic field lines and changes in the magnetic 
field topology. While these phenomena are behind the emission 
of large amounts of energy, they are strictly forbidden within the 
ideal-MHD limit due to magnetic flux conservation and so can not 
be studied employing this limit. 

Besides having considerably smaller conductivities, low- 
density higly magnetized plasmas are present rather generically 
around magnetized objects, constituting what is referred to as 
the "magnetosphere". In such regions magnetic stresses are much 
larger than magnetic pressure gradients and cannot be properly bal- 
anced; as a result, the magnetic fields have to adjust themselves 
so that the magnetic stresses vanish identically. This scenario is 
known as the force free regime (because the Lorentz force vanishes 
in this case) and while the equations governing i t can be seen as the 
low-inertia limit o f the ideal-MHD equations l lKomissarovl 1200 j ; 
lMcKinney||2006bh . the force-free limit is really distinct from the 



1 Note that this expression for the electrical conductivity is roughly correct 
for densities in the range 10 10 — 10 14 g cm -3 and temperatures in the range 
10 6 — 10 8 K, but provide s a reasonable estimat e also at larger temperatures 
of ~ 10 9 - 10 10 K [cf. . |PotekhinetaIUl999h l. 



ideal-MHD one. This represents a considerable complication since 
it implies that it is usually not possible to decribe, within the same 
set of equations, both the interior of compact objects and their mag- 
netospheres. 

Theoretical work to derive a fully relativistic theory of non- 
ideal hydrodynamics and non-ideal MHD has been carried out by 
several authors in the past I Israelii 19761 : IStewartll 19771 ; ICarteJl99ll ; 
lLichnerowicdll96llAniielll989h and is particularly simple in the 
case of the resistive MHD description. The purpose of this work 
is indeed that of proposing the solution of the relativistic resistive 
MHD equations as an important step towards a more realistic mod- 
elling of astrophysical plasmas. There are a number of advantages 
behind such a choice. First, it allows one to use a single mathe- 
matical framework to describe both regions where the conductiv- 
ity is large (as in the interior of compact objects) and small (as in 
magnetospheres), and even the vacuum regions outside the compact 
objects where the MHD equations trivially reduce to the Maxwell 
equations. Second, it makes it possible to account self-consistently 
for those resistive effects, such as current sheets, which are energet- 
ically important and could provide a substantial modification of the 
whole dynamics. Last but not least, the numerical solution of the 
resistive MHD equations provides the only way to control and dis- 
tinguish the physical resistivity from the numerical one. The latter, 
which is inevitably present and proportional to truncation error, is 
also completely dependent on the specific details of the numerical 
algorithm employed and on the resolution used for the solution. 

As noted already by several authors, the numerical solution 
of the ideal-MHD equations is considerably less challenging than 
that of the resistive MHD equations. In this latter case, in fact, 
the equations become mixed hyperbolic-parabolic in Newtonian 
physics or hyperbolic with stiff relaxation terms in special rela- 
tivity. The presence of stiff terms is the natural consequence of the 
fact that the diffusive effects take place on timescales that are in- 
trinsically larger than the dynamical one. Stated differently, in such 
equations the relaxation terms can dominate over the purely hyper- 
bolic ones, posing severe constraints on the timestep for the evolu- 
tion. While considerable work has already been made to introduce 
numeric al techniques to ach i eve efficient implementa ti ons in either 

20071 ; iKomissarovl 



regim e i Komissarovl |2004|: 



120071 ; iRevnolds et alj|200a 



Komissarov et al 



Graves et al. 2008[), the use of these 



techniques in fully three-dimensional simulations is still difficult 
and expensive. 

In order to benefit from the many advantages discussed above 
in the use of the resistive MHD equations, we here present a novel 
approach for the solution of the relativistic resistive MHD equa- 
tions exploiting the properties of implicit-explicit (IMEX) Runge 
Kutta methods. This approach represents a simple but effective so- 
lution to the problem of the vastly different timescales without sac- 
rificing either simplicity in the implementation or the numerical 
efficiency. By examining a number of tests we illustrate the ac- 
curacy of our approach under a variety of conditions and demon- 
strate its robustness. In additi on, we also c ompa re it with the al- 
ternative method proposed by IKomissarovl d2007l) for the solution 
of the same set of relativistic resistive MHD equations. This latter 
approach employs Strang-splitting techniques and the analytical in- 
tegration of a reduced form of Ampere's law. While it works well 
in a number of cases, it has revealed to be unstable when applied 
to discontinuous flows with large conductivities; such difficulties 
were not encountered when solving the same problem within the 
IMEX implementation. 

Because our approach effectively treats within a unified 
framework both those regions of the flow which are fluid-pressure 
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dominated and those which are instead magnetic-pressure domi- 
nated, it could find a number of applications and serve as a first 
step towards a more realistic modeling of relativistic astrophysical 
plasmas. 

Our work is organized as follows. In Sect. [2] we present the 
system of equations describing a resistive magnetized fluid, while 
in Section[3]we discuss the problems related to the numerical evo- 
lution of this system of equations and the numerical approaches 
developed to solve them. In particular, we introduce the basic fea- 
tures of the IMEX Runge-Kutta schemes and recall their stability 
properties. In Sect [4] we instead explain in detail the implementa- 
tion of the IMEX scheme to the resistive MHD equations. Finally, 
in Sect.[5]we present the numerical tests carried out either in one or 
two dimensions and that span several prescriptions for the conduc- 
tivity. Section|5]is also dedicated to the comparison with the Strang- 
splitting technique. The conclusions and the perspectives for future 
improvements are presented in Sect. [6] while Appendix lAl reviews 
our space discretization of the equations. 

Hereafter we will adopt Gaussian units such that c = 1 and 
employ the summation convention on repeated indices. Roman in- 
dices a, b, c, ... are used to denote spacetime components (i.e., from 
to 3), while i,j, k, ... are used to denote spatial ones; lastly, bold 
italics letters represent vectors, while bold letters represent tensors. 



Levi-Civita symbol e abc is non-zero only for spatial indices. Note 
that the electromagnetic fields have no components parallel to n a 
(i.e., E a n a = = B a n a ). 

By using the decomposition of the Maxwell tensor l[3j, the 
equations ([T}-(|2j can be split into directions which are parallel and 
orthogonal to n a to yield the familiar Maxwell equations 



V-E = q, (4) 

V-.B = 0, (5) 

d t E-X7 x B = J, (6) 

d t B + V x E = , (7) 



where we have decomposed also the current vector I a = qn a + J a , 
with q being the charge density, qn a the convective current and J a 
the conduction current satisfying J a n a = 0. 

The current conservation equation d a I a = follows from the 
antisymmetry of the Maxwell tensor and provides the evolution of 
the charge density q 

d t q + V • J = , (8) 

which can be obtained also directly by taking the divergence of ® 
when the constraints {4]l-l[5} are satisfied. 



2 THE RESISTIVE MHD DESCRIPTION 

An effective description of a fluid in the presence of electromag- 
netic fields can be made by considering three different sets of equa- 
tions governing respectively the electromagnetic fields, the fluid 
variables and the coupling between the two. In particular, the elec- 
tromagnetic part can be described via the Maxwell equations, while 
the conservation of energy and momentum can be used to express 
the evolution of the fluid variables. Finally, Ohm's law, whose exact 
form depends on the microscopic properties of the fluid, expresses 
the coupling between the electromagnetic fields and the fluid vari- 
ables. In what follows we review these three sets of equations sepa- 
rately, discuss how they then lead to the resistive MHD description, 
and how the latter reduces to the well-known limits of ideal-MHD 
and of the Maxwell equations in vacuum. Our presentation will be 
focussed on the special-relativistic regime, but the extension to gen- 
eral relativity is rather straightforward and will be presented else- 
where. 



The Maxwell equations 

Th e special relativistic M axwell equations can be written 
as dLandau & LifshitJ 198Ch 



d b F ab = r 

d b *F ab = , 



(1) 

(2) 

where F ab and *F ab are the Maxwell and the Faraday tensor re- 
spectively and I a is the electric current 4- vector. A highly-ionized 
plasma has essentially zero electric and magnetic susceptibilities 
and the Faraday tensor is then simply the dual of the Maxwell ten- 
sor. This tensor provides information about the electric and mag- 
netic fields measured by an observer moving along any timelike 
vector n a , namely 



b 771a . abc j-t 

n E + e B c 



(3) 



We are considering n a to be the time-like traslational killing vector 
field in a flat (Minkowski) spacetime, so n a = (—1, 0, 0, 0) and the 



The hydrodynamic equations 

The evolution of the matter follows from the conservation of the 
stress-energy tensor 



dtT = , 

and the conservation of baryon number 

d a (pu a ) = , 



(9) 



(10) 



where p is the rest-mass density (as measured in the rest frame of 
the fluid) and u a is the fluid 4-velocity. The stress-energy tensor 



T ab describing a perfect fluid minimally coupled to an electromag- 

(11) 



netic field is given by the superposition 

-ifluici . mem 



rri rT-itlUlCl . rrv 

J-ab — J-ab ~T ± a b > 



where 



rj-iO 



F ac F b -i(F cd F cd )g ab , 



T<ab i a b i ab 

fluid = hu u +pg . 



(12) 
(13) 



Here h = p(l + e) + p is the enthalpy, with p the pressure and e 
the specific internal energy. 

The conservation law (|9j can be split into directions parallel 
and orthogonal to n a to yield the familiar energy and momentum 
conservation laws 



OtT + V ■ F T = , 
d t S + V • F s = , 



(14) 
(15) 



where we have introduced the conserved quantities {r, S}, which 
are essentially the energy density r = T a t,n a n b and the energy flux 



\{E 2 + B 2 ) + hW 2 -p, 



density Si = T a in a , and whose expressions are given by 

(16) 

S = E x B + hW z v . (17) 

Here v is the velocity measured by the inertial observer and W = 
—n a u a = — v 2 is the Lorentz factor. The fluxes can then be 
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written as 
F T = E x B + hW 2 V , 
F s ee - SB + hW 2 vv + 



\{E 2 + B 2 ) 



(18) 
(19) 



Finally, the conservation of the baryon number J 1 0| > reduces to the 
continuity equation written as 

d t D + V ■ F D = , (20) 

where we have introduced another conserved quantity D = pW 
and its flux F d ee pWv. 



Ohm's law 

As mentioned above, Maxwell equations are coupled to the fluid 
ones by means of the current 4- vector I a , whose explicit form will 
depend in general on the electromagnetic fields and on the local 
fluid properties. A standard prescription is to consider the current 
to be proportional to the Lorentz force acting on a charged particle 
and the electrical resistivity r\ to be a scalar function. Ohm's law, 
written in a Lorentz invariant way, then reads 



la + (I b U b )u a = CT Fab U ■ 



(21) 



with g ee I/77 being the electrical conductivity of the medium. 
Expressing \2\\ in terms of the electric and magnetic fields one 
obtains the familiar form of Ohm's law in a general inertial frame 



J = a W[E + v x B - (E ■ v)v] + q v 



(22) 



Note that the conservation of the electric charge <[8j provides the 
evolution equation for the charge density q (i.e., the projection of 
the 4-current J along the direction n), while Ohm's law provides a 
prescription for the (spatial) conduction current J (i.e., the compo- 
nents of I orthogonal to n). 

It is important to recall that in deriving expression l !22t for 
Ohm's law we are implicitly assuming that the collision frequency 
of the constituent particles of our fluid is much larger that the 
typical oscillation frequency of the plasma. Stated differently, the 
timescale for the electrons and ions to come into equilibrium is 
much shorter than any other timescale in the problem, so that no 
charge separation is possible and the fluid is globally neutral. This 
assumption is a key aspect of the MHD approximation. 

The well-known ideal-MHD limit of Ohm's law can be ob- 
tained by requiring the current to be finite even in the limit of in- 
finite conductivity (a — ► 00). In this limit Ohm's law i22l then 
reduces to 



E + v x B - (E ■ v)v = . 



(23) 



Projecting this equation along v one finds that the electric field does 
not have a component along that direction and then from the rest of 
the equation one recovers the well-known ideal-MHD condition 



E = —v x B . 



(24) 



stating that in this limit the electric field is orthogonal to both B and 
v. Such a condition also expresses the fact that in ideal MHD the 
electric field is not an independent variable since it can be be com- 
puted via a simple algebraic relation from the velocity and magnetic 
vector fields. 

Summarizing: the system of equations of the relativistic re- 
sistive MHD approximation is given by the constraint equations 
(|4}-(|5]l, evolution equations (|6j-([8]l, dT4b — <TT3T) and l [20l . where the 
fluxes are given by Eqs. l|18|l-l|19t and the 3-current is given by 



Ohm's law i22\ . These equations, together with a equation of state 
(EOS) for the fluid and a reasonable model for the conductivity, 
completely describe the system under consideration provided con- 
sistent initial and boundary data are defined. 



Different limits of the resistive MHD description 

At this point it is useful to point out some properties of the rela- 
tivistic resistive MHD equations discussed so far, to underline their 
purely hyperbolic character and to contrast them with those of other 
forms of the resistive MHD equations which contain a parabolic 
part instead. To do this within a simple example, we adopt the New- 
tonian limit of Ohm's law ( 122b . 



J = a[E + v X B) 



(25) 



where we have neglected terms of order 0(v 2 /(?), obtaining the 
following potentially stiff equation for the electric field 



d t E - V x B = -a\E + v x B] 



(26) 



Assuming now a uniform conductivity and taking a time derivative 
of Eq. Q, we obtain the following hyperbolic equation with re- 
laxation terms (henceforth referred simply as hyperbolic-relaxation 
equation) for the magnetic field 



1 



[d u B - V 2 B] = [d t B - V x (v x B)} . 



(27) 



If the displacement current can be neglected, i.e., dtE ~ 
dttB ~ 0, equation J27t reduces to the familiar parabolic equa- 
tion for the magnetic field 



d t B - V x (v x B) V 2 B = . 

a 



(28) 



where the last term is responsible for the diffusion of the magnetic 
field. It is important to stress the significant difference in the charac- 
teristic structure between equations d27t and J28l >, Both equations 
reduce to the same advection equation in the ideal-MHD limit of 
infinite conductivity (a — > 00) indicating the flux-freezing condi- 
tion. However, in the opposite limit of infinite resistivity (a — > 0) 
Eq. l |28t tends to the (physically incorrect) elliptic Laplace equa- 
tion V 2 .B = while Eq. l |27b reduces to the (physically correct) 
hyperbolic wave equation for the magnetic field. 



2.1 The augmented MHD system 

The set of Maxwell equations described above can also be cast in 
an extended fashion which includes two additional fields, ip and 
<j>, introduced to control dynamically the constraints of the system, 
i.e., Eqs I0 and (O. This "augmented" system reads 

d b {F ab + ipg ab ) = r-ra/m 11 , (29) 
db( *F ab + 4>g ab ) = -K<f)n a . (30) 

Clearly, the standard Maxwell equations ((TJ-lO are recovered 
when tp = <f> — and we are in this way extending the space of 
solutions of the original Maxwell equations to include those with 
non-vanishing {ip, </>}. 

The evolution of these extra scalar fields can be obtained by 
taking a partial derivative d a of the augmented Maxwell equations 
<|29b — J30b and using the antisymmetry of the Maxwell and Faraday 
tensors together with the conservation of charge to obtain 



d a d a i> 
d a d a <t> 



-Kd a (ipn a ) , 
-n,d a ((f>n a ) . 



(3D 
(32) 
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It is evident that these represent wave equations with sources for 
the scalar fields {ip, <f>}, which propagate at the speed of light while 
being damped if k > 0. In particular, for any positive n, they de- 
cay exponentially over a timescale ~ 1/k to the trivial solution 
ip = 4> = and the augmented system then reduces to the stan- 
dard Maxwell equations, including the constraints © and l|5}. This 
approach, n amed hyperbolic div ergence cleaning in the context of 
ideal MHD jPedner et alj|2002h . was proposed as a simple way of 
solving the Maxwell equations and enforcing the conservation of 
the divergence-free condition for the magnetic field. 

Ado pting this approach and following the formulation pro- 
posed by IKomissarovl J2007I) . the evolution equations of the aug- 
mented Maxwell equations fl29M30ll can then be written as 



dt^ + V-E = q-rnjj, (33) 

d t cj> + V-B = - K <t>, (34) 

d t E - V x B + Vtp = —J , (35) 

d t B + V x E + V4> = . (36) 



The system of equations d33M 36ll. together with the current con- 
servation l[8j, is the one we will use for the numerical evolution 
of the electromagnetic fields within the set of relativistic resistive 
MHD equations. 



3 EVOLUTION OF HYPERBOLIC-RELAXATION 
EQUATIONS 

While the ideal-MHD equations are well suited to an efficient nu- 
merical implementation, the general system of relativistic resistive 
MHD equations brings about a delicate issue when the conductivity 
in the plasma undergoes very large spatial variations. In the regions 
with high conductivity, in fact, the system will evolve on timescales 
which are very different from those in the low-conductivity region. 
Mathematically, therefore, the problem can be regarded as a hyper- 
bolic one with stiff relaxation terms which requires special care to 
capture the dynamics in a stable and accurate manner. In the next 
Section we discuss a simple example of a hyperbolic equation with 
relaxation which exhibits the problems discussed above and then 
introduce implicit-explicit (IMEX) Runge Kutta methods to deal 
with these kind of equations. In essence, these methods treat the 
advection character of the system with strong-stability preserving 
(SSP) explicit schemes, while the relaxation character with an L- 
stable diagonally implicit Runge Kutta (DIRK) scheme. After pre- 
senting the scheme, its properties and some examples, we discuss 
in detail its application to the resistive MHD equations. 

3.1 Hyperbolic systems with relaxation terms 

A prototypical hyperbolic equation with relaxation is given by 

d t U = F{U) + -R(U) , (37) 

where e > is the relaxation time (not necessarily constant either 
in space or in time), F(U) gives rise to a quasilinear system of 
equations (i.e., F(U) depends linearly on first derivatives of U), 
and R does not contain derivatives of U . 

In the limit e — > oo (corresponding for the resistive MHD 
equations to the case of vanishing conductivity) the system is hy- 
perbolic with propagation speeds bounded by c/,. This maximum 
bound, together with the length scale L of the system, define a 
characteristic timescale t>j = L/ch of the hyperbolic part. In the 



opposite limit e —> (corresponding to the case of infinite conduc- 
tivity), the system is instead said to be stiff, since the timescale e of 
the relaxation (or stiff) term R(U) is in general much larger than 
the timescale th of the hyperbolic part F(U). In such a limit, the 
stability of an explicit scheme is only achieved with a timestep 
size At ^ e. This requirement is certainly more restrictive than the 
Courant-Lewy-Friedrichs (CFL) stability condition At ^ Ax/ch 
for the hyperbolic part and makes an explicit integration impracti- 
cal. The development of efficient numerical schemes for such sys- 
tems is challenging, since in many applications the relaxation time 
can vary by several orders of magnitude across the computational 
domain and, more importantly, to much beyond the one determined 
by the speed Ch- 

When faced with this issue several strategies can be adopted. 
The most straightforward one is to consider only the stiff limit e — » 
0, where the system is well approximated by a suita ble reduced 
set of conservation laws called "equilibrium system" dChen et all 
1 1994 such that 

R(U) = 0, (38) 
dtU = G(U). (39) 

where U is a reduced set of variables. This approach can be fol- 
lowed if the resulting system is also hyperbolic. This is precisely 
the case in the resistive MHD equations for vanishing resistivity 
rj — ► (or (T — > oo). In this case, the equations reduce to those of 
ideal MHD and describe indeed an "equilibrium system" in which 
the magnetic field is simply adverted with the flow. As discussed 
earlier, this limit is often adequate to describe the behaviour of 
dense astrophysical plasmas, but it may also stray away in the mag- 
netospheres. A more general approach could consist of dividing 
the computational domain in regions in each of which a simplified 
set of equations can be adopted. As an example, the ideal-MHD 
equations could be solved in the interior of compact objects, the 
force-free MHD equations could be solved in the magnetosphere, 
and finally the Maxwell equations for the vacuum regions outside 
the compact object. However, this approach requires the overall 
scheme to suitably match the different regions so as to obtain a 
global solution. This task, unfortunately, is far from being straight- 
forward and, to date, it lacks a rigorous definition. 

An alternative approach consists of considering the origi- 
nal hyperbolic-relaxation system in the whole computational do- 
main and then employ suitable numerical schemes that work 
for all reg ions. Among such schemes is the Strang-splitting 
tec hnique dStrand 1 1968b . which has been recently applied 
by IKomissarovl ( 120071) for the solution of the (special) relativis- 
tic resistive MHD equations. The Strang-splitting scheme provides 
second-order accuracy if each step is at least second-order accu- 
rate, and this property is maintained under su itable assumptions 
even for stiff problems jjahnke & Lubichll2000l) . In practice, how- 
ever, higher-order accuracy is difficult to obtain even in non-stiff 
regimes with this kind of splitting. Moreover, when applied to hy- 
perbolic systems with relaxation, Strang-splitting schemes reduce 
to first-order accuracy since the kernel of the relaxation operator is 
non-trivial and corresponds to a singular matrix in the linear case, 
therefore invalidating the assumptions made by I Jahnke & Lubichl 

2 Implicit schemes could avoid this issue at an increased computational 
cost; however, an explicit second order accurate method approaching iter- 
atively the Crank-Nicholson scheme has been shown, in a simple model 
with hyperbolic-relaxation terms, to work well when dealing with smooth 
profiles without being too costly (M. Choptuik, private communication) 
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to ensure high-order accuracy. iKomissarovl J2007h avoided 
this problem by solving analytically the stiff part in a reduced 
form of Ampere's law. Although this procedure works well for 
smooth solutions, our implementation of the method has revealed 
problems when evolving discontinuous flows (shocks) for large- 
conductivities plasmas. Moreover, it is unclear whether the same 
procedure can be adopted in more general configurations, where an 
analytical solution may not be available. 

As an alternative approach to the methods solving the rela- 
tivistic resistive MHD equations on a single computati onal domain, 
we here in t roduce an IMEX Runge-Kutta method jAsher et"al] 
Il995l 1 19971 : |Pareschill200ll : IPareschi & Russol fcoOS ) to cope with 
the stiffness problems discussed above. These methods, which are 
easily implemented, are still under development and have few (rel- 
atively minor) drawbacks. The most serious one is a degradation to 
first or second-order accuracy for a range of values of the relaxation 
time e. However, since High-Resolution Shock-Capturing (HRSC) 
schemes usually employed for the solution of the hydrodynamic 
equations already suffer from similar effects at discontinuities, the 
possible degradation of the IMEX schemes does not spoil the over- 
all quality numerical solution when employed in conjunction with 
HRSC schemes. The next sections review in some detail the IMEX 
schemes and our specific implementation for the relativistic resis- 
tive MHD equations. 



3.2 The IMEX Runge-Kutta methods 

The IMEX Runge-Kutta schemes rely on the application of an im- 
plicit discretization scheme to the stiff terms and of an explicit one 
to the non-stiff ones. Wh en applied to system d37t it takes the form 
jPareschi & Russoll2005l) 

U {i) = U n + AtJ2 a t3 F(U u) )+AtJ2 a l3 -R(U u) ) , 
3=1 3=1 e 

U n+1 = U n + AtJ2^F(U (l) )+AtJ2^-R(U { 



(40) 

where U^ 1 ' are the auxiliary intermediate values of the Runge- 
Kutta scheme. The matrices A = (dij) and A — (dij) are 
v x v matrices such that the resulting scheme is explicit in F 
(i.e., hij = for j i) and implicit in R. An IMEX Runge-Kutta 
scheme is characterized by these two matrices and the coefficient 
vectors Wj and Wj. Since simplicity and efficiency in solving the 
implicit part at each step is important, it is natural to consider di- 
agonally implicit Runge-Kutta (DIRK) schemes (i.e., cnj = for 
j > i) for the stiff terms. 

A particularly convenient way of describing an IMEX Runge- 
Kutta scheme is offered by the Butcher notation, in which the 
scheme is by a double tableau of the type jButchedl 19871 . [2003h 



(41) 



c 


A 


c 


A 






T 

UJ 



where the index T indicates a transpose and where the coefficients 
c and c used for the treatment of non-autonomous systems are given 
by 



i-l 

Ci = ^ ] dij , 
3=1 



= E 

3=1 



(42) 



Table 1. Tableau for the explicit (left) implicit (right) IMEX-SSP2 (2, 2, 2) 
L-stable scheme 












7 


7 





1 


1 





1-7 


I-27 


7 




1/2 


1/2 




1/2 


1/2 



7=1- 



1 



restrictions in some of the coefficients of their respective Butcher 
tableaus. Although each of them separately can have an arbitrary 
accuracy, this does not ensure that the combination of the two 
schemes will preserve the same accuracy. In addition to the above 
conditions for each Runge-Kutta scheme, there are also some addi- 
tional conditions combining terms in the two tableaus which must 
be fulfilled in order to achieve a global accuracy order for the com- 
plete IMEX scheme. 

Since the details of these methods are not widely known, we 
first consider a simple example to fix ideas. A second-order IMEX 
scheme can be written in the tableau form given in TableQ] The in- 
termediate and final steps of this IMEX Runge-Kutta scheme would 
then be written explicitly as 



jjW = Jjr. 

u {2) = xr 



AtF(U 
At, 



in. 



(2)-, 



u 



The accuracy of each of the Runge-Kutta is achieved by imposing 



+ — [(1 - 2^)R(U W ) + yR(U 
e 

"+ 1 = [/" + ^[^(f/W) + F(U (2) )] 

+ ^[R(U W ) + R(U^)]. 

Note that at each sub-step an implicit equation for the auxiliary 
intermediate values U' 1 ' must be solved. The complexity of invert- 
ing this equation will clearly depend on the particular form of the 
operator R(U). 



3.2.1 Stability properties of the IMEX schemes 

Stable solutions of conservation-type equations are usually ana- 
lyzed in terms of a suitable norm being bounded in time. With U n 
representing the solution vector at the time t = n At, then a se- 
quence {U n } is said to be "strongly stable" in a given norm || • | 
provided that \\U n+1 \\ «C ||f7 n || for all n ^ 0. 

The most commonly used norms for analyzing schemes for 
nonlinear systems are the Total- Variation (TV) norm and the in- 
finity norm. A numerical scheme that maintains strong stability at 
the discrete leve l is ca lled Strong Stability Preserving (SSP) (see 
ISpiteri & Ruuthl j2002l) for a detailed description of optimal SSP 
schemes and their pro perties). Because of the s tability properties 
of the IMEX schemes l lPareschi & Russo|[20o3) . it follows that if 
the explicit part of the IMEX scheme is SSP, then the method is 
SSP for the equilibrium system in the stiff limit. This property is 
essential to avoid spurious oscillations during the evolution of non- 
smooth data. 

The stability of the implicit part of the IMEX scheme is en- 
sured by requiring that the Runge-Kutta is "L-stable" and this rep- 
resents an essential condition for stiff problems. In practice, this 
amounts to requiring that the numerical approximation is bounded 
in cases when the exact solution is bounded. A more strict defini- 
tion can be derived starting from a linear scalar ordinary differential 
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Table 2. Tableaux for the explicit (first row) and implicit (second row) 
IMEX SSP-schemes. We use the standard notation SSPfc(s, cr, p), where 
k denotes the order of the SSP scheme and the triplet (a, a, p) characterizes 
respectively the number of stages of the implicit scheme (s), the number of 
stages of the explicit scheme (cr), and the order of the IMEX scheme (p). 



SSP2 (3,3,2) 














u 


("I 
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1/2 


1/2 
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1/3 
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1/3 
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SSP3 (3,3,2) 
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1 


1 





n 
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1/2 


1/4 


1/4 


n 
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SSP3 (4,3,3) 
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1/6 


2/3 


a 


a 













—a 
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1 





1 — OL 
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1/2 
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V 


1/2 - - 







1/6 




1/6 







a 



a = 0.24169426078821 , /3 ■ 
7 = l-l/%/2, v z - 



2/3 

0.06042356519705 , 
0.12915286960590. 



equation, namely 



(43) 



In this case it is easy to define the stability (or amplification) func- 
tion C(z) as the ratio of the solutions at subsequent timesteps 
C{z) = * n+1 /* n , where z = Atq. A Runge-Kutta scheme is 
then said to be L-stable if |C(z )| < 1 (i.e., it is bounded) and 
C(oo) = teutcheJl98ll2003h . 

There are a number of IMEX Runge-Kutta schemes available 
in the literature and we report here only some of the second and 
third-order schemes which satisfy the condition that in the limit 
e — > 0, the solution correspond s to that of the equilibrium sys- 
tem l l38l ( IPareschi & Russol2005h . These are giv en in their Butcher 
tablea u form in Table [2] and are taken from IPareschi & Russol 
d2005l) . In all these schemes the implicit tableau corresponds to 
an L-stable scheme. The tableaus are reported in the notation 
SSPfc(s, a, p), where k denotes the order of the SSP scheme and 
the triplet (s, o,p) characterizes respectively the number of stages 
of the implicit scheme (s), the number of stages of the explicit 
scheme (<r), and the order of the IMEX scheme (p). 



4 IMEX RUNGE-KUTTA SCHEME FOR THE 
AUGMENTED RESISTIVE MHD EQUATIONS 

Having reviewed the main properties of the IMEX schemes, we 
now apply them to the particular case of the special relativistic re- 
sistive MHD equations. Our goal is to consider a numerical imple- 
mentation of the general system that can deal with standard hydro- 
dynamic issues (like shocks and discontinuities) as well as those 
brought up by the stiff terms discussed in the previous Section. 
Hence, we adopt high-resolution shock-capturing algorithms (see 
AppendixfA} together with IMEX schemes. Because the first ones 
involve the introduction of conserved variables in order to cast the 
equations in a conservative form, we first discuss how to implement 
the IMEX scheme within our target system and subsequently how 
to perform the transformation from the conserved variables to the 
primitive ones. 



4.1 IMEX schemes for the Maxwell-Hydrodynamic 
equations and treatment of the implicit stiff part 

For our target system of equations it is possible to introduce a 
natural decomposition of variables in terms of those whose evo- 
lution do not involve stiff terms and those which do. More specif- 
ically, with the electrical resistivity r\ playing the role of the re- 
laxation parameter e, the vector of fields U can be split in two 
subsets {X,Y}, with X = {E} containing the stiff terms, and 
Y = {B, tp, </>, q, r, S, D} the non-stiff ones. 

Following the prototypical Eq. {37J, the evolution equations 
for the relativistic resistive MHD equations can then be schemati- 
cally written as 



d t Y = 
d t X = 



F Y (X,Y), 
F X (X,Y) + 



e(Y) 



R X (X,Y) 



(44) 
(45) 



where the relaxation parameter e is allowed to depend also on the 
Y non-stiff fields. The vector Y can be evolved straightforwardly 
as it involves no stiff term. We further note that for our particular 
set of equations, it is convenient to write the stiff part as 



Bx(X,Y) = A(Y)X + S x (Y) 



(46) 



As a result, the procedure to compute each stage U^' of the IMEX 
scheme can be performed in two steps: 

(i) Compute the explicit intermediate values { X * , Y * } from all 
the previously known levels, that is 



(47) 



X* = X n + AtY l *» F x(U^) 

3=1 



.At£igL^(E/«>), 
(48) 

where we have defined e^' = e(~F^) and atj /e"' in Eq. {48} is 
a simple division and not a contraction on dummy indices. 

(ii) Compute the implicit part, which involves only X, by solv- 
ing 



= Y* , 

= X * +At^rR x (U ( ' ) ) 



(49) 
(50) 



Note that the implicit equation, with the previous assumption l |46| >. 
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can be inverted explicitly 

x d) = M (Y*)(X* +au ^L S X (Y*)), (51) 

M(Y*) = [I — an ^^A(Y*)]~ 1 , (52) 

since the form of the matrix [/ — an AtA(Y*) / 'e^'] is known 
explicitly in terms of the evolved fields. 

The explicit expressions for stiff part are then given simply by 

R E = -WE + W (E ■ v)v -Wv x B , (53) 
S E = -Wv x B , (54) 

with the matrix A defined as 

(-l+vl V X Vy V X V z \ 

v x v y -1 + Vy VyV z . (55) 

V z V x V z V y -1 + v\ j 

Hence, the matrix M can be computed explicitly to obtain 

1 / a + W + aW 2 vl aW 2 v x v y aW 2 v x v z \ 

— I aW 2 v x v y a + W + aW 2 v 2 aW 2 v y v z 
m \ aW 2 v z v x aW 2 v z v v a + W + aW 2 v 2 J 

where m = W 2 a + Wa 2 + W + a and a = an a (i) At. 

Summarizing: First, an intermediate state {E*} is found 
through the evolution of the non-stiff part for the electric field. Sec- 
ond, if the velocity v is known, the evolution of the stiff part can be 
performed by acting with M to obtain 

E = M(v) [E* + an At a (l) S E {v, B)] . (56) 

At this point the approach proceeds with the conversion from the 
conserved variables to the primitive ones. Because of the coupling 
between the electric and the velocity fields, such a procedure is 
rather involved and more complex than in the ideal-MHD case; a 
detailed discussion of how to do this in practice will be presented 
in Sect. [42] 

It is interesting to highlight the consistency at two known lim- 
its of the implicit solution of the stiff part. In the ideal-MHD limit 
(i.e., g — > oo) the first term of Eq. l |56t vanishes, while the contri- 
bution of the second term leads to the ideal-MHD condition J24t . 
On the other hand, in the vanishing conductivity limit (i.e., a — > 0) 
the second term in Eq. d56| l vanishes, and the matrix reduces to the 
identity one M(v) — I. In this case, the electric field is obtained 
only by evolving the explicit part, i.e., E = E* . 

Finally, it is important to stress that one could, in principle, 
have considered the alternative route of adopting instead X = 
{E, q}, so that the right-hand-side of q would be considered stiff 
with R q — and S q — V • R E . However, this choice could 
lead to spurious numerical oscillations in the solution since the 
fluxes of q can be discontinuous, while they would be evolved 
with an implicit Runge-Kutta. As it has been shown under fairly 
gener al conditions, high-o rder SSP schemes are necessarily ex- 
plicit dGottlieb et al j|200lh . so it follows that this part of the equa- 
tions cannot be evolved with the implicit Runge-Kutta unless a low- 
order scheme is implemented. 



4.2 Transformation of conserved variables to primitive ones 

As mentioned in the previous Section, in order to evolve our sys- 
tem of equations, the fluxes {F T ,Fs , F d} must be computed 
at each timestep. These fluxes depend on the primitive fields 



{p, p, v, E, B}, which must be recovered from the evolved 
conserved fields {D, t, S, E, B}. These quantities are related 
by complicated equations which become transcendental except for 
particularly simple equations of state (EOS). As a result, the con- 
version must be in general pursued numerically and the primitive 
variables are then given by the roots of the function 



f(p) = P{p, e) - : 



(57) 



where p(p, e) is given by the chosen EOS and p is the trial value 
for the pressure eventually leading to the primitive variables. 

Note that since = Y* [cf. Eq. l|49]l], the values of the 
conserved quantities {D, r, S, B} at time (n + l)At are ob- 
tained by evolving their non-stiff evolution equations which, how- 
ever, provide only an approximate solution for the electric field 
{E*}. As discussed in the previous Section, the final solution for 
the electric field E requires the inversion of an implicit equation 
and, hence, is a function of the velocity v and of the fields {B , E* } 
[cf. Eq. J56H . However, the velocity is a primitive quantity and thus 
not known at the time (n+ 1) At. It is clear, therefore, that it is nec- 
essary to obtain, at the same time, the evolution of the stiff part of 
the equations and the conversion of the conserved quantities into to 
the primitive ones. In what follows we describe how to do this in 
practice using an iterative procedure. 

(i) Adopt as initial guess for the velocity its value at the previous 
time level v — v n . The electric field E is computed by Eq. l |56t as 
a function of [E* , v,B). 

(ii) Adopt as initial guess for the pressure its value at the previ- 
ous time level p — p n . Compute in the following order 

S E x B 



W 



r-{E 2 + B 2 )/2+p 
1 



D 

P = W' 

t-(E 2 + B 2 )/2 -DW+p(l- W 2 ) 
DW 



(58) 



(iii) Solve numerically Eq. J57t by means of an iterative 
Newton-Raphson solver, so that the solution at the iteration m + 1 
can be computed as 

/(Pra) 



Om+l = Pm 



f'(Pf 



(59) 



The derivative of the function f(p) needed for the Newton- 
Raphson solver can be computed as 

f(p) = v 2 c 2 ~-l, (60) 

with c 3 being the local speed of the fluid which, for an ideal-fluid 

EOS p(p, e) = (r — 1) p e is given by 

(iv) With the newly obtained values for the velocity v and the 
pressure p, the steps (i)-(iii) can be iterated until the difference 
between two successive values falls below a specified tolerance. 

The approach discussed above is a simple procedure that can 
be implemented straightforwardly and works well for moderate 
ratios of \B\ 2 /p, converging in less than 10 iterations both for 
smooth electromagnetic fields and for discontinuous ones. Faster 
and more robust procedures to obtain the primitive variables cer- 
tainly ca be implemented, but this is beyond the scope of this work. 
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Figure 1. Magnetic field component B y for a large-amplitude CP Alfven 
wave and for three different resolutions Ax = {1/50, 1/100, 1/200}. 
The conductivity is constant with a magnitude of <r = 10 6 . The agreement 
betweem the exact solution and that corresponding to the high resolution 
one is excellent. 




-1 



X 

Figure 2. Magnetic field component B y in a self-similar current sheet. The 
solution is computed with N = 200 gridpoints (Aa; = 1/200) and is 
shown at the initial time t = 1 and at t = 10. The conductivity is uniform 
with a magnitude of a = 10 2 (i.e., r\ = 1/cr = 0.01). The numerical 
solution is in excellent agreement with the exact one. 




5 NUMERICAL TESTS 

In this section we present several one-dimensional (ID) or two- 
dimensional (2D) tests which have been used to validate the im- 
plementation of the IMEX Runge-Kutta schemes in the different 
regimes of relativistic resistive MHD. In all these tests we employ 
the ideal-fluid EOS with V = 2 for the ID tests and T = 4/3 in 
the 2D ones. The different tests span several prescriptions for the 
conductivity and compare the solutions obtained either with those 
expected in the ideal-MHD limit or with those obtained with the 
Strang-splitting technique. 

More specifically, in ID we consider large-amplitude circu- 
larly polarized (CP) Alfven waves to test the ability of the code to 
reproduce the ideal-MHD results when adopting a very large con- 
ductivity. The intermediate conductivity regime is instead tested by 
simulating a self-similar current sheet. Finally, a large range of uni- 
form and non-uniform conductivities are used for a representative 
shock-tube problem. In 2D, on the other hand, we first consider 
a commonly employed test for ideal-MHD codes corresponding 
to a cylindrical explosion. Subsequently, we simulate a toy model 
for a "magnetized neutron star" when modelled as a cylindrically 
symmetric density distribution obeying a Gaussian-profile. The be- 
haviour of the magnetic field is studied again for a range of constant 
and non-uniform conductivities. 



5.1 One-dimensional tests 

5.1.1 Large amplitude CP Alfven waves 

This test is discussed in detail by iDel Zanna et ail d2007l) and we 
report here only a short summary. The solution describes the propa- 
gation of a large amplitude circularly-polarized Alfven waves along 
a uniform background field Bo in a domain with periodic boundary 
conditions. The exact solution in the i deal-MHD limit and a ssum- 
ing v x = for simplicity, is given by jDel Zanna et al 1l2007h 



where B x = Bo, k is the wave vector, r)A is the amplitude of the 
wave and the special relativistic Alfven speed va is given by 



(By, B z ) = tiaBq (cos[k(x - v A t)],sin[k(x - v A t)\) 
(v y ,v z ) = (B y ,B z ) , 



2 

va 



h + BUl + vl) 



1r)ABl 



h + Bl{i + n l) 



(62) 



(63) 

In practice, using such ideal-MHD solution it is possible to assess 
the accuracy of evolution of the resistive equations by requiring 
that for very large conductivities the numerical solution approaches 
the exact one as the resolution is progressively increased. It is also 
worth remarking that although we do not expect the solution of 
the resistive MHD equations to converge to that of ideal MHD for 
any finite value of a, we also expect the differences between the 
two to be 0(v/a) and thus negligibly small for sufficiently large 
values. For this reason, we have performed the evolution with a 
high uniform conductivity of a — 10 6 for three different reso- 
lutions N = {50, 100, 200} covering the computational domain 
x 6 [—0.5, 0.5]. In addition, the initial data parameters have been 
chosen so that p = p — tja = 1 and Bo = 1.1547, thus yielding 
Va = 1 /2, with a full period being achieved at t = 2. 

Fig. [TJ confirms this expectation by reporting the component 
By after one period and thus overlapping with the initial one (at 
t = 0) for the highest resolution. This test shows clearly that in the 
limit of very high conductivity the resistive MHD equations tend 
to a solution which is very close to the same solution obtained in 
the ideal-MHD limit. The convergence rate measured for the differ- 
ent fields is consistent with the second-order spatial discretization 
being used as expected for smooth flows (see Appendix [At. 



5.7.2 Self-similar current sheet 

The details of this test are described by iKomissarovl < |2007l) . so 
again we provide here only a short description for completeness. 
We assume that the magnetic pressure is much smaller than the 
fluid pressure everywhere, with a magnetic field given by B — 
(0, B y (x, t),Q), where B y (x, t) changes sign within a thin current 
layer of width Al. Provided the initial solution is in equilibrium 
(p = const.), the evolution is a slow diffusive expansion of the 
layer due to the resistivity and described by the diffusion equation 
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Figure 3. Left panel: Magnetic field component B y in the solution of the shock-tube problem. Different lines refer to three different resolutions and to the 
exact ideal-MHD solution at t = 0.4. The conductivity is uniform with a magnitude of <ro = 10 6 . Right panel: The same as in the left panel but for different 
uniform conductivities. Note that for ctq = the solution describes a discontinutiy propagating at the speed of light and corresponding to Maxwell equations 
in vacuum. As the conductivity increases, the solution tends to the ideal-MHD one. 
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Figure 4. Differences in the magnetic field component B y between the 
numerical solution computed with either the Strang or the IMEX schemes 
and the exact solution of the shock-tube in the ideal-MHD limit. The differ- 
ences are computed for several uniform conductivities, although the Strang- 
splitting technique does not yield a stable solution for values larger than 
cro ~ 7000 for the reference resolution of Ax = 1/400 (i.e. with 400 
gridpoints). Shown in the inset is the maximum conductivity for which a 
solution was possible, <r max , as a function of the number of gridpoints, N. 

[cf. Eq. {28]! with v = 0] 

d t B y - -d 2 B y = . (64) 
a 

As the system expands, the width of the layer becomes much larger 

than AZ and it evolves in a self-similar fashion. For t > 0, the 
analytical exact solution is given by 

B y (x,t) = B erf Q^f) > ( 65 > 

where £ = t/x 2 and "erf" is the error function. This solution 
ca n be used for testi ng the moderate resistive regime. Follow- 
ing lKomissarovl d2007t) . and in order to avoid the singular behaviour 
at t = 0, we have chosen as initial data the solution at t = 1 with 
p = 50, p = 1, E = v = and a = 100. The domain covers the 
region x £ [—1.5, 1.5] with N = 200 points. 



The numerical simulation is evolved up to t = 10 and then 
the numerical and the exact solution are compared in Fig. [2] The 
two solutions match so well that they are not distinguishable on 
the plot, thus, showing that the intermediate-conductivity regime is 
also well described by our method. 

5.1.3 Shock-tube problem 

As prototypical shock-tub e test we consider a simple MHD version 
of the Brio and Wu test terio & wj|l988r) , where the initial left 
and right states are separated at £ = 0.5 and are given by 

(p L ,p L ,B y L ) = (1.0,1.0,0.5), 
(p R ,p R ,B R ) = (0.125,0.1,-0.5), 

while all the other fields set to 0. We consider both uniform and 
non-uniform conductivities. In the latter case we adopt the follow- 
ing prescription 

a = a D~< , (66) 

thus allowing for nonlinearities in the dependence of the conductiv- 
ity on the conserved quantity D. This is one of the simplest cases, 
but in realistic situations a more general expression for the conduc- 
tivity can be assumed, where a is a function of both the rest-mass 
density and of the specific internal energy, i.e., o — a(p, e). 

The exac t solution of the ideal MHD Riemann problem 
was found by jGiacomazzo & Re zzolla 2006]), and in our partic- 
ular case_jt_Jias_beeri_comput a publicly available code 
Isee lGiacomazzo &Rezzollal ( f2006|) l. When B x = 0, the structure 
of the solution contains only two fast waves, a rarefaction mov- 
ing to the left and a shock moving to the right, with a tangential 
discontinuity between them. More demanding Riemann problems 
have also been performed but the procedure to convert the con- 
served variables into the primitive ones has shown in these case a 
lack of robustness for large ratios of \B\ 2 /p. 

We have first considered the case of uniform (7 = 0) and 
very large conductivity (<jo = 10 6 ) as in this case we can use the 
solution in the ideal-MHD limit as a useful guide. The profile of 
the magnetic field component B y for three different resolutions 
Ax = {1/100, 1/200, 1/400} and the exact solution are shown 
in the left panel of Fig. [3] at t = 0.4. Overall, the results indicate 
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Figure 5. Left panel: Evolution of a non-uniform conductivity a in the shock-tube problem for different values of 7 and indicated by the different lines 
(o"0 = 10 6 for all lines). Notice the large variability on the magnitude of the conductivity. Right panel: The same as in the left panel but for the magnetic field 
component B y . 



that even in the presence of shocks our numerical solution of the 
resistive MHD tends to the ideal-MHD solution as the resolution 
is increased. It is also interesting to study the behaviour of the so- 
lution for different values of the constant ao while still keeping a 
uniform conductivity (i.e., 7 = 0). This is shown in the right panel 
of FiglU which displays the different solutions obtained, and where 
it is possible to see how they change smoothly from a wave-like so- 
lution for cr = to the ideal-MHD one for cr = 10 6 . 

This set up is also useful to perform a comparison between the 
IMEX and the Strang-splitting approaches. In Fig. [4] we show the 
L 1 -norm of the difference between the numerical solution obtained 
with both schemes and the ideal-MHD exact solution, for different 
values of the conductivity with TV = 400 points. 

Several comments are in order. Firstly, the reported difference 
between the numerical solution for the resistive MHD equations 
and the ideal-MHD equations should not be interpreted as an er- 
ror given that the latter is not the correct solution of the equa- 
tions. Hence, the fact that the use of a Strang-splitting method 
yields smaller differences is simply a measure of its ability of better 
capture steep gradients. Secondly, while the IMEX approach does 
not show any sign of instability for 00 ranging between 10 2 and 
10 9 , the implementation adopting the Strang-splitting technique 
becomes unstable for moderately high values of the conductivity 
and, at least for the shock-tube problem, no numerical solution was 
possible for ao > 7000 at the above resolution. Increasing the reso- 
lution can help increase the maximum value of the resistivity which 
can be handled, but since this gain is only linear with the number 
of gridpoints aiming for higher conductivities results impractical. 
This is shown in the inset of Fig. [4] which reports the maximum 
conductivity for which a solution was possible, <r max , as a function 
of the number of gridpoints, N. Finally, we note that the difference 
between the IMEX numerical solution and the exact ideal-MHD 
one saturates between ao ~ 10 5 — 10 6 . This is not surprising since 
the differences are expected to be 0(1/ a), and thus the saturation 
in the differences essentially provides a measure of our truncation 
error at the resolution used. 

A more challenging test is offered by the solution of the shock- 
tube in the presence of a non-uniform conductivity. In particular, we 
have considered the same initial states and the same non-uniform 
conductivity discussed above, but used different values for the ex- 
ponent 7 in l |66t while keeping ao constant. The results of this test 



are shown in the left panel of Fig. [5] where the conductivity is plot- 
ted at t — 0.4 for several values of 7. Note that the conductivity 
traces the evolution of the rest-mass density and that the solution 
can be found also when a varies of almost 12 orders of magni- 
tude across the grid. Similarly, the right panel of Fig.[5]displays the 
component B y for the different values of 7. It should be stressed 
that because of the relation i66l between a and p, the region on 
the left has at this time a very high conductivity and the numerical 
solution tends to the ideal-MHD one. The opposite happens on the 
right region, where the conductivity is lower for higher values of 7. 
Clearly, the results presented in Fig. [5] show that our implementa- 
tion can handle non-uniform (and quite steep) conductivity profiles 
even in the presence of shocks. 

5.2 Two-dimensional tests 

5.2.1 The cylindrical explosion 

We now consider problems involving shocks in more than one di- 
mension. A demanding test for the relativistic codes is the cylin- 
drical blast wave expanding in a plasma with an initially uniform 
magnetic field. Although there is no exact solution for this prob- 
lem, strong symmetric explosions are useful tests since shocks are 
present in all the possible directions and the numerical implementa- 
tion is therefore tested in all of its parts. For this test we set a square 
domain (x,y) € [— 6, 6] with a resolution Aa; = Ay = 1/200. 
The initial data is such that inside the radius r < 0.8 the pressure 
is set to p — 1 while the density to p — 0.01. In the intermedi- 
ate region 0.8 ^ r ^ 1.0 the two quantities decrease exponen- 
tially up to the exterior region r > 1, where the ambient fluid has 
p = p = 0.001. The magnetic field is uniform with only one non- 
trivial component B = (0.05, 0, 0). The other fields are set to be 
zero (i.e., E — q = 0), which is consistent within the ideal-MHD 
approximation. 

The evolution is performed with a high conductivity a = 10 6 
in order to recover the solution from the ideal-MHD approxima- 
tion. As shown in Fig. [6] which reports the magnetic field compo- 
nents B x (left panel) and B y (right panel) at time t — 4, we obtain 
results tha t are qualitatively similar to those published in differ- 
ent works jKomissarovll999al : iNeilsen et alj200ol iDel Zanna et all 
|2007| J Komissaro vl l2007l) . While a strict comparison with an exact 
solution is not possible in this case, the solution found matches ex- 



© 2008 RAS, MNRAS 000. [Wl5\ 



12 Palenzuela, Lehner, Reula, Rezzolla 




4.77e-03 2.82s>01 



Figure 6. Magnetic field components B x (left panel) and B y 

tremely well the one obtained with another 2D code solving the 
ideal MHD equations. Most importantly, however, the figure shows 
that the solution is regular everywhere and that similar results can 
be obtained also with smaller values of the conductivity (e.g., no 
significant difference was seen forcr > 10 4 ). 

5.2.2 The cylindrical star 

We next consider a toy model for a star, thought as an infinite 
column of fluid aligned with the z-axis but with compact support 
in other directions. Because of the symmetry in the z-direction, 
d,U — for all the fields and the problem is therefore two- 
dimensional. More specifically, we consider initial data given by 

p = Po e- (r/r ° )2 , (67) 
v = (v r ,v*,v z )=p(0,u>*,0), (68) 

B = (B r ,B*,B*)=p (o,0,2B„(l-£)\ , (69) 

where r = \Jx 2 + y 2 is the cylindrical radial coordinate. The 
other fields can be computed at the initial time by using the poly- 
tropic EOS p = p r , the ideal-MHD expression ( 124b for the elec- 
tric field, and the electric charge from the constraint equation 
q = V ■ E. We have chosen r a = 0.7, p = 1.0, ^ = 0.1 and 
Bo = 0.05. An atmosphere ambient fluid with p — 0.01 is added 
outside the cylinder. Finally, the resolution is Ax = 1 /200 and the 
domain is (x, y) G [—3, 3]. 

This simple problem exhibits some of the issues present in a 
magnetized rotating neutron star: a compactly supported rest-mass 
density distribution, an azimuthal velocity field and a poloidal mag- 
netic field. Suitable source terms describing a gravitational poten- 
tial have been added to the Euler equations in order to get, at least 
at the initial time, a stationary solution. In the ideal-MHD limit the 
magnetic lines are frozen in the fluid and thus a static profile is also 
expected for the magnetic field. 

In the left panel of Fig. [7] we plot the slice y — of the mag- 
netic field component B z at t = 14 as obtained from the evolution 
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(right panel) for the cylindrical explosion test at time t = 4. 

of the resistive MHD system for different uniform conductivities in 
the range ao e [10 2 , 10 6 ]. In the limiting case ao = the solution 
corresponds to a wave propagating at the speed of light (i.e., the 
solution of the Maxwell equations in vacuum), while for large val- 
ues of ao the solution is stationary (as expected in the ideal-MHD 
limit). The behaviour observed in the left panel Fig. [7] is also the 
expected one: the higher the conductivity, the closer the solution 
is to the stationary solution of the ideal-MHD limit. For low con- 
ductivities, on the other hand, there is a significant diffusion of the 
solution, which is quite rapid for ao < 10 2 and for this reason those 
values are not plotted here. We note that values of the conductivity 
larger than ao > 10 7 lead to numerical instabilities that we believe 
are coming from inaccuracies in the evolution of the charge density 
q, and which contains spatial derivatives of the current vector. In 
addition, the stiff quantity E x is seen to converge only to an order 
~ 1.5. This can be due to the "final layer" problem of the IMEX 
methods, which is known to produce a degradation on the accuracy 
of the stiff quantities. Luckily, this does not spoil the convergence 
of the non-stiff fields, which are instead second-order convergent. 
It is possible that the use of stiffly-accurate schemes can solve this 
degradation of the convergence and this is an issue we are presently 
exploring. 



We finally consider the same test, but now employing the non- 
uniform conductivity given by Eq. l |66t with ao = 10 6 and dif- 
ferent values for 7. The results are presented in the right panel of 
Fig. [7] which shows that the magnetic fields inside the star are ba- 
sically the same in all the cases, stressing the fact that the interior 
of the star will not be significantly affected by the exterior solution, 
which has much smaller conductivity. However, the electromag- 
netic fields outside the star do change significantly for different 
values of 7, underlining the importance of a proper treatment of 
the resistive effects in those regions of the plasma where the ideal- 
MHD approximation is not a good one. 
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6 CONCLUSIONS 

We have introduced Implicit-Explicit Runge-Kutta schemes to 
solve numerically the (special) relativistic resistive MHD equations 
and thus deal, in an effective and robust way, with the problems in- 
herent to the evolution of stiff hyperbolic equations with relaxation 
terms. Since for these methods the only limitation on the size of the 
timestep is set by the standard CFL condition, the approach sug- 
gested here allows to solve the full system of resistive MHD equa- 
tions efficiently without resorting to the commonly adopted limit of 
the ideal-MHD approximation. 

More specifically, we have shown that it is possible to split 
the system of relativistic resistive MHD equations into a set of 
equations that involves only non-stiff terms, which can be evolved 
straightforwardly, and a set involving stiff terms, which can also 
be solved explicitly because of the simple form of the stiff terms. 
Overall, the only major difficulty we have encountered in solving 
the resistive MHD equations with IMEX methods arises in the con- 
version from the conserved variables to the primitive ones. In this 
case, in fact, there is an extra difficulty given by the fact that there 
are four primitive fields which are unknown and have to be inverted 
simultaneously. We have solved this problem by using extra itera- 
tions in our ID Newton-Raphson solver, but a multidimensional 
solver is necessary for a more robust and efficient implementation 
of the inversion process. 

With this numerical implementation we have carried out a 
number of numerical tests aimed at assessing the robustness and 
accuracy of the approach, also when compared to other equiva- 
len ts ones, such as the Strang-splitting method recently proposed 
by iKomissarovl J2007h . All of the tests performed have shown 
the effectiveness of our approach in solving the relativistic resis- 
tive MHD equations in situations involving both small and large 
uniform conductivities, as well as conductivities that are allowed 
to vary nonlinearly across the plasma. Furthermore, when com- 
pared with the Strang-splitting technique, the IMEX approach has 
not shown any of the instability problems that affect the Strang- 
splitting approach for flows with discontinuities and large conduc- 
tivities. 

While the results presented here open promising perspectives 
for the implementation of IMEX schemes in the modelling of rel- 
ativistic compact objects, at least two further improvements can be 
made with minor efforts. The first one consists of the generalization 



of the (special) relativistic resistive MHD equations with a scalar 
isotropic Ohm's law to the general relativistic case, and its appli- 
cation to com pact astrophysical bodies such a mag netized binary 
neutron stars l lAnderson et al J[2008l ; iLiu et alj|2008l) . The solution 
of the resistive MHD equations can yield different results not only 
in the dynamics of the magnetosphere produced after the merger, 
but also provide the possibility to predict, at least in some approx- 
imation, the electromagnetic radiation produced by the merger of 
these objects. The second improvement consists of considering a 
non-scalar and anisotropic Ohm's law, so that the behaviour of the 
currents in the magnetosphere can be described by using a very 
high conductivity along the mag netic lines and a ne gligibly small 
one in the transverse directions (Komissarov 2004). Such an im- 
provement may serve as a first step towards an alternative mod- 
elling of force-free plasmas. 
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APPENDIX A: TVD SPACE DISCRETIZATION 

We are generically interested in solving hyperbolic conservation 
laws of the form 

d t U + d k k F{U) = S(U) , (Al) 

where U is the vector of the evolved fields, k F are their fluxes 
and S contains the sources terms. The semi-discrete version of this 
equation, in one dimension, is simply given by 

QtUi = _ F I+1/2 -F^ 1/2 + sm ^ 
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where ^±1/2 are consistent numerical fluxes evaluated at the in- 
terfaces between numerical cells. These consistent fluxes are com- 
puted by using HRSC methods, which are based on the use of Rie- 
mann solvers. More specifically, we have implemented a modifica- 
tion of the Local Lax-Friedr ichs approximate Riemann solver in- 
troduced bv lAlic etal.1d2007t) , which only needs the spectral radius 
(i.e., the maximum eigenvalue) of the system. In highly relativistic 
cases, like the ones we are interested in, the spectral radius is close 
to the light speed c = 1 and so the Local Lax-Friedrichs reduces to 
the simpler Lax-Friedrichs flux 



F, 



i+1/2 



= t:[Fl+F r + (ul - Ur) 



(A3) 



where ul,ur are the reconstructed solutions on the left and on 
the right of the interface and Fl,Fr their corresponding fluxes. 
The standard procedure is then to reconstruct the solution ul, ur 
by interpolating with a polynomial and then compute the fluxes 
Fl = F(ul) and Fr — F(ur). In ou r implementation we first 
recombine the fluxes and the solution as dAlic et al.ll2007l) 



Fi ±Ui 



(A4) 



Then, using a piecewise linear reconstruction, these combinations 
can be computed on the left/right of the interface as 



Fr — F i+1 - - A i+1 , 



(A5) 



where A* are just the slopes used to extrapolate F^ to the inter- 
faces. Finally, the consistent flux is computed by a simple average 

1, 



Fi+i 



/2 



-\ F L 



Fc 



For a linear reconstruction the slopes can be written as 

A+ = L(F+ 1 -F+,F+-F+ 1 ), 



L(F^_ 2 F i , 1 ,F i , 1 F^~) , 



(A6) 



(A7) 



so that it is trivial to check that the standard Lax-Friedrichs ( I A3b is 
recovered when A+ = A~ . The choice of these slopes becomes 
crucial in the presence of shocks or very sharp profiles, while the 
use of some nonlinear operators L(x, y) preserves the Total Varia- 
tion Diminishing (TVD) condition on the interpolating polynomial. 
In this way, the TVD schemes capture accurately the dynamics of 
strong shocks without the oscillations which appear with standard 
finite-difference discretizations. Monotonicity is typically enforced 
by making use of slope limiters and we have in particular imple- 
mented the Monotonized Centered (MC) limiter 

L(x,y) = i[sign(s:)+sign(j/)]min(2|x|,2|y|, , (A8) 

which provides a good compromise between robustness and accu- 
racy. Note that with linear reconstruction the scheme is second- 
order accurate in the smooth regions, although it drops to first order 
near shocks and at local extrema. 
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